This document presents preliminary analyses characterizing the NDVI of certain places and census tracts throughout the Denver area.
We obtained the normalized difference vegetation index (NDVI) from the Landsat-8 satellite at a resolution of 30 meters squared at a roughly 15-day interval between 2016 and 2021 for the following area around Denver:
We gathered this data using the rgee package, an R package that facilitates connection to the Google Earth Engine via Python. Please see these two scripts for specific details on the process for gathering and processing the NDVI data in this area:
https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1a_get_landsat_ndvi_denver.R
Here is the NDVI over the study area on a cloud-free day, July 4, 2021:
#Update 5/26/22 this may not work with mapview anymore, so use plot instead for speed and
setwd(here("data-processed"))
pal_terrain_col = rev(terrain.colors(100))
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif")
ndvi_den_metro_terr_5_yr$`20210704_NDVI` %>% plot()

#old code:
# mv_good_date =ndvi_den_metro_terr_5_yr$`20210704_NDVI` %>%
# raster::raster() %>%
# mapview(
# col.regions = pal_terrain_col,
# layer.name = "NDVI, July 4, 2021")
# mv_good_date
First, pick a few places where we would expect NDVI to be high in the summer on a cloud-free day. If it’s not high, then we can assume the image is bad (i.e., clouds in the way). These test places were determined by examining the cloud-free day (July 4, 2021) above. We chose three areas: an area east of Evergreen, Colorado; a plot in City Park; and a plot in the Indian Tree Golf Club:
(These areas are defined here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R)
setwd(here("data-processed"))
getwd()
## [1] "/Users/michaeldgarber/Dropbox/Work/CSU-UCSD post-doc/post-doc-research/green-space-denver-project/green-space-denver/data-processed"
load("bbox_evergreen_east.RData")
load("bbox_indian_tree_golf.RData")
load("bbox_city_park.RData")
load("den_metro_bbox_custom.RData")
mv_evergreen_east = bbox_evergreen_east %>%
mapview(col.regions = "red", layer.name = "Evergreen East")
mv_indian_tree_golf = bbox_indian_tree_golf %>%
mapview(col.regions = "red", layer.name = "Indian Tree Golf")
mv_city_park =bbox_city_park %>%
mapview(col.regions = "red", layer.name = "City Park plot")
mv_den_metro_bbox_custom = den_metro_bbox_custom %>%
mapview(layer.name = "Study area", alpha.regions = .2)
mv_den_metro_bbox_custom+ mv_evergreen_east + mv_indian_tree_golf + mv_city_park
Let’s look at the temporal NDVI trend for these three places, without excluding any dates. This is the mean NDVI on that day. That is, each place includes several pixels of size 30 meters squared. This is the average of the NDVI for those pixels for that place on that day.
setwd(here("data-processed"))
load("ndvi_test_places_day_wrangle.RData")
ndvi_test_places_day_wrangle %>%
ggplot(
aes(
x=date,
y=ndvi_mean
))+
geom_line(aes(colour = test_place_name))+
geom_point(aes(colour = test_place_name))+
scale_x_date(labels=date_format("%Y-%m-%d"), date_breaks = "3 months")+
theme(axis.text.x = element_text(angle = 45, hjust=1))

Hmm, if we look closely, we can see that there are some very low NDVI values for some of these places even in the summer, which is not plausible. Those low values must indicate obstruction by clouds. For example, look at a day when NDVI was measured as very low in City Park despite it being mid-May when we would expect it to be higher. So I’m going to exclude dates like these.
setwd(here("data-processed"))
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif")
ndvi_den_metro_terr_5_yr$`20180517_NDVI` %>% #Use 2018-05-17
plot()

#old code that uses mapview. decide to use plot instead. see abobve
# mv_bad_date =ndvi_den_metro_terr_5_yr$`20180517_NDVI` %>% #Use 2018-05-17
# raster::raster() %>%
# mapview(
# col.regions = pal_viridis_trunc,
# layer.name = "NDVI, May 17, 2018")
# mv_bad_date
The specific NDVI thresholds for each thresholds are noted in this script (https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R) and warrant further discussion.
The valid dates are restricted to May, June, July, and August and specifically are:
setwd(here("data-processed"))
load("lookup_date_is_valid_all.RData") #load dataset of valid dates created elsewhere
lookup_date_is_valid_all %>%
filter(date_is_valid_all==1) %>%
dplyr::select(date) %>%
pull()
## [1] "2016-05-08" "2016-06-01" "2016-06-09" "2016-06-17" "2016-07-03"
## [6] "2016-08-12" "2017-05-09" "2017-06-02" "2017-06-26" "2017-07-12"
## [11] "2017-07-28" "2017-08-13" "2017-08-21" "2017-08-29" "2018-05-25"
## [16] "2018-06-10" "2018-06-18" "2018-06-26" "2018-07-04" "2018-07-12"
## [21] "2018-08-05" "2018-08-13" "2018-08-21" "2018-08-29" "2019-05-09"
## [26] "2019-06-26" "2019-07-04" "2019-08-21" "2019-08-29" "2020-05-24"
## [31] "2020-06-25" "2020-07-27" "2020-08-12" "2020-08-28" "2021-05-25"
## [36] "2021-06-02" "2021-06-10" "2021-07-04" "2021-07-28" "2021-08-29"
Load the places-of-interest dataset, which is created here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_read_denver_native_zones.R
setwd(here("data-processed"))
load("places_native_geo.RData")
load("places_native_nogeo.RData")
load("native_places_ndvi_day_nogeo.RData")
Print information about the places of interest, including approximate percent nativity where available.
places_native_nogeo %>%
dplyr::select(place_name, native_percent, place_type, place_id) %>%
kable()
| place_name | native_percent | place_type | place_id |
|---|---|---|---|
| Denver Botanic Gardens, 100% Native | 100 | native spectrum | 1 |
| Plains Conservation Center, 100% Native | 100 | native spectrum | 2 |
| Green Mountain Park, 85% Native | 85 | native spectrum | 3 |
| Hogback along C-470, 75% Native | 75 | native spectrum | 4 |
| Suburban Open Space 1, Chatfield H.S., 50% Native | 50 | native spectrum | 5 |
| Suburban Open Space 2, Columbine Hills Church, 30% Native | 30 | native spectrum | 6 |
| Denver Botanic Gardens at Chatfield, 10% Native | 10 | native spectrum | 7 |
| City Park, 0% Native | 0 | native spectrum | 8 |
| Chatfield Meadow Restoration | NA | high diversity | 9 |
| Chatfield Prairie Garden | NA | high diversity | 10 |
| City Park Greenhouses | NA | high diversity | 11 |
| Denver Botanic Gardens Green Roof | NA | high diversity | 12 |
| Kenderick Lake Xeriscape Garden | NA | high diversity | 13 |
Map these places by type (nativity spectrum or high-plant diversity).
During valid dates only (above).
The first set of polygons sent over.
native_places_ndvi_day_nogeo %>%
filter(date_is_valid_all==1) %>%
filter(place_type == "native spectrum") %>%
ggplot(aes(
x=date,
y=ndvi_med # median
))+
geom_ribbon( #ribbon around 25th and 75th percentile
aes(ymin =ndvi_25, ymax = ndvi_75 ),
alpha=.4
)+
ylab("NDVI, Median") +
xlab("Date") +
geom_line( size=.7 ) +#note size better than lwd
geom_point()+
scale_x_date(breaks = pretty_breaks())+
scale_y_continuous(
limits = c(0, NA), #force axis origin to be zero
breaks= seq(0,0.8,by = 0.1))+
theme(
axis.text.x=element_text(angle=60, hjust=1),
panel.border = element_rect(
colour = "gray72", size = 0.5, fill=NA))+
facet_grid( # facet them
cols = vars(place_name_fac),
labeller = label_wrap_gen() #wrap facet labels
)

This is the second set of polygons sent.
native_places_ndvi_day_nogeo %>%
filter(date_is_valid_all==1) %>%
filter(place_type == "high diversity") %>%
ggplot(aes(
x=date,
y=ndvi_med # median
))+
geom_ribbon(
aes(ymin =ndvi_25, ymax = ndvi_75 ),
alpha=.4
)+
ylab("NDVI, Median") +
xlab("Date") +
geom_line( size=.7 ) +
geom_point()+
scale_x_date(breaks = pretty_breaks())+
scale_y_continuous(
limits = c(0, NA),
breaks= seq(0,0.8,by = 0.1))+
theme(
axis.text.x=element_text(angle=60, hjust=1),
panel.border = element_rect(
colour = "gray72", size = 0.5, fill=NA))+
facet_grid(
cols = vars(place_name_fac),
labeller = label_wrap_gen() #wrap facet labels
)

Among the first set of polygons with values for percent nativity.
native_places_ndvi_day_nogeo %>%
filter(date_is_valid_all==1) %>%
filter(place_type == "native spectrum") %>%
ggplot(aes(
x=native_percent,
y=ndvi_med # median
))+
geom_point(
aes(colour=place_name_fac), size = 1.5,
alpha=.5 #varying alpha to illustrate density
) +
xlab("Percent native") +
ylab("NDVI, median") +
scale_y_continuous(
limits = c(0, NA), #force axis origin to be zero
breaks= seq(0,0.8,by = 0.1))+
scale_color_hue(
name = "Place name"
)

setwd(here("data-processed"))
load("places_native_nogeo.RData")
load("native_places_ndvi_day_nogeo.RData")
native_places_ndvi_day_nogeo %>%
filter(date_is_valid_all==1) %>%
group_by(place_type, place_name_fac) %>%
summarise(
ndvi_mean = mean(ndvi_mean_wt, na.rm=TRUE), # mean of means
ndvi_25th = quantile(ndvi_med, probs =c(0.25), na.rm=TRUE), #percentile of medians
ndvi_med = median(ndvi_med, na.rm=TRUE), #percentile of medians
ndvi_75th = quantile(ndvi_med, probs =c(0.75), na.rm=TRUE) ) %>% #percentile of medians
ungroup() %>%
left_join(places_native_nogeo, by = c("place_name_fac", "place_type")) %>%
dplyr::select(
starts_with("place_type"),
starts_with("place_name_fac"),
starts_with("area_ft2"),
starts_with("area_mi"),
starts_with("ndvi_mea"),
starts_with("ndvi_25"),
starts_with("ndvi_med"),
starts_with("ndvi_75th"),
) %>%
as_tibble() %>%
knitr::kable(
booktabs = TRUE,
col.names = c("Place type", "Place name", "Area (sq. ft.)", "Area (sq. mi.)", "NDVI, mean", "NDVI, 25th-ile", "NDVI, 50th-ile", "NDVI, 75th-ile"),
digits = c(0, 0, 0, 3, 2, 2, 2, 2, 2)
)
| Place type | Place name | Area (sq. ft.) | Area (sq. mi.) | NDVI, mean | NDVI, 25th-ile | NDVI, 50th-ile | NDVI, 75th-ile |
|---|---|---|---|---|---|---|---|
| high diversity | Chatfield Meadow Restoration | 4391 | 0.000 | 0.48 | 0.47 | 0.50 | 0.50 |
| high diversity | Chatfield Prairie Garden | 2922 | 0.000 | 0.37 | 0.35 | 0.37 | 0.37 |
| high diversity | City Park Greenhouses | 9596 | 0.000 | 0.31 | 0.27 | 0.28 | 0.28 |
| high diversity | Denver Botanic Gardens Green Roof | 12882 | 0.000 | 0.34 | 0.34 | 0.36 | 0.36 |
| high diversity | Kenderick Lake Xeriscape Garden | 35418 | 0.001 | 0.40 | 0.36 | 0.37 | 0.37 |
| native spectrum | Denver Botanic Gardens, 100% Native | 21585 | 0.001 | 0.54 | 0.52 | 0.55 | 0.55 |
| native spectrum | Plains Conservation Center, 100% Native | 793905 | 0.028 | 0.16 | 0.23 | 0.26 | 0.26 |
| native spectrum | Green Mountain Park, 85% Native | 83536011 | 2.996 | 0.42 | 0.36 | 0.40 | 0.40 |
| native spectrum | Hogback along C-470, 75% Native | 1293880 | 0.046 | 0.36 | 0.32 | 0.35 | 0.35 |
| native spectrum | Suburban Open Space 1, Chatfield H.S., 50% Native | 161894 | 0.006 | 0.40 | 0.37 | 0.41 | 0.41 |
| native spectrum | Suburban Open Space 2, Columbine Hills Church, 30% Native | 198596 | 0.007 | 0.36 | 0.26 | 0.31 | 0.31 |
| native spectrum | Denver Botanic Gardens at Chatfield, 10% Native | 97135 | 0.003 | 0.53 | 0.40 | 0.54 | 0.54 |
| native spectrum | City Park, 0% Native | 2664873 | 0.096 | 0.67 | 0.67 | 0.70 | 0.70 |
In this section, we summarize NDVI for census tracts in Denver County. Before characterizing NDVI of each census tract, we removed bodies of water. NDVI values of water approach -1, and we did not want to “penalize” a census tract for having large amounts of water.
Here is a map of the median NDVI of the census tracts in the area on 2021-07-04:
Here is a histogram of the distribution of median NDVI for these census tracts on June 10, 2021:
den_co_tract_ndvi_day_geo %>%
filter(date == "2021-06-10") %>%
ggplot(aes(ndvi_med))+
geom_histogram()

In this section, we examine the NDVI levels of public green space in this area, including city parks and nature reserves. Please see this code for more details on how we gathered the polygons for the greens spaces: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_load_denver_osm_parks_water.R
Here are the parks we downloaded:
setwd(here("data-processed"))
load("den_jeff_co_green_space_public.RData")
den_jeff_co_green_space_public %>%
mapview(zcol = "osm_value", layer.name = "Green space type")
Like with census tracts, we removed bodies of water from these green spaces before measuring NDVI. Here is the NDVI on June 10, 2021:
setwd(here("data-processed"))
load("den_metro_green_space_ndvi_day_geo.RData")
load("den_jeff_co_geo.RData")
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
mv_den_jeff_co_geo = den_jeff_co_geo %>%
mapview(
layer.name = "County name",
color = c("red", "orange"),
col.regions = c("red", "orange"),
lwd=2,
zcol = "county_name_short", alpha.regions = 0)
mv_den_metro_green_space_ndvi_day_geo = den_metro_green_space_ndvi_day_geo %>%
dplyr::select( #limit some of the variables before mapview
starts_with("osm"),
starts_with("ndvi"),
starts_with("date"),
starts_with("county")) %>%
filter(date == "2021-06-10") %>%
mapview(
layer.name = "NDVI, Median",
col.regions = pal_viridis_trunc,
zcol = "ndvi_med"
)
mv_den_metro_green_space_ndvi_day_geo+mv_den_jeff_co_geo
Most parks are in the city, and most protected areas are in the mountains. Summarize NDVI by type of green space:
setwd(here("data-processed"))
library(knitr)
load("den_metro_green_space_ndvi_day_nogeo.RData")
den_metro_green_space_ndvi_day_nogeo %>%
filter(date_is_valid_all==1) %>%
group_by(osm_value) %>%
summarise(
ndvi_25th = quantile(ndvi_med, probs =c(0.25), na.rm=TRUE),
ndvi_med = median(ndvi_med, na.rm=TRUE),
ndvi_75th = quantile(ndvi_med, probs =c(0.75), na.rm=TRUE),
) %>%
ungroup() %>%
knitr::kable()
| osm_value | ndvi_25th | ndvi_med | ndvi_75th |
|---|---|---|---|
| nature_reserve | 0.3313458 | 0.4128266 | 0.4128266 |
| park | 0.3722610 | 0.4735859 | 0.4735859 |
| protected_area | 0.3428352 | 0.4256270 | 0.4256270 |
Copyright © 2022 Michael D. Garber